
This is the documentation of the MODULE \texttt{LapackAPI}, a set
of routines to perform calls to the Linear Algebra PACKage. 



\section{Subroutine \texttt{GEVP(A, B,L [,V,Is])}}
\index{GEVP@Subroutine \texttt{GEVP(A, B,L [,V,Is])}}

\subsection{Description}

The routine \texttt{GEVP}, solves the generalized eigenvalue problem 
\begin{displaymath}
  Av^{(n)} = \lambda_{(n)}B v^{(n)}
\end{displaymath}

\subsection{Arguments}

\begin{description}
\item[\texttt{A, B}:]  Real or complex double precision two
  dimensional array. The matrices of the general eigenvalue problem.

\item[\texttt{L(:): }] Real double precision one dimensional array. A
  vector containing the left generalized eigenvalues.

\item[\texttt{U(:,:):}] Real or complex double precision two dimensional
  array. Output. Optional. The left eigenvectors as columns.

\item[\texttt{Is}] Integer. Output. Optional. 0 if the computation was
  sucessfull, other integer in any other case. 
\end{description}

\subsection{Example}

\begin{lstlisting}[emph=GEVP,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Using lapack library to solve a generalized
                   eigenvalue problem.,
                   label=gevp]
Program GEVPprg

  USE NumTypes
  USE Linear
  USE LapackAPI

  Integer :: N = 2, Is
  Complex (kind=DPC), Allocatable :: A(:,:), B(:,:), v(:,:)
  Real (kind=DP), Allocatable :: l(:), R(:)

  ALLOCATE(A(N,N), B(N,N), v(N,N), l(N), R(N))

  A(1,1) = 1.0
  A(1,2) = 2.0
  A(2,1) = 2.0
  A(2,2) = 3.0

  B(1,1) = 2.0
  B(1,2) = 3.0
  B(2,1) = 3.0
  B(2,2) = 8.0

  Write(*,*)'Det(B)', Det(B)

  CALL GEVP(A, B, l, v, Is)
  Write(*,*)'Success: ', Is
  Write(*,*)l
  Write(*,*)
  Do I = 1, 2
     Do J = 1, 2
        Write(*,*)I,J, v(I,J)
     End Do
  End Do
  Write(*,*)

  Write(*,*)'Residuals: '
  Do J = 1, N
     Do I = 1, N
        R(I) = Abs(Sum( (A(I,:)-l(J)*B(I,:))*v(:,J)))
     End Do
     Write(*,*)R
  End Do


  Stop
End Program GEVPprg
\end{lstlisting}

\section{Subroutine \texttt{SVD(A, L [, U, VT])}}
\index{SVD@Subroutine \texttt{SVD(A [, U, VT]) }}

\subsection{Description}

The routine \texttt{SVD}, performs the singular value decomposition of
the matrix $A(:,:)$, optionally returning the right and left
eigenvectors. 

\subsection{Arguments}

\begin{description}
\item[\texttt{A}:]  Real double precision two dimensional array. The
  matrix whose singular value decomposition we want to know.

\item[\texttt{L(:): }] Real double precision one dimensional array. A
  vector containing the singular values of $A$.

\item[\texttt{U(:,:):}] Real double precision two dimensional
  array. Output. Optional. The left eigenvectors as columns.

\item[\texttt{VT(:,:):}] Real double precision two dimensional
  array. Output. Optional. The right eigenvectors as columns.
\end{description}

\subsection{Example}

\begin{lstlisting}[emph=SVD,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Using lapack library to perform a SVD decomposition.,
                   label=svd]
Program SVDtest

  USE NumTypes
  USE LapackAPI

  Real (kind=DP), Allocatable :: A(:,:), U(:,:), VT(:,:), Sv(:),&
       & SA(:,:)

  Real (kind=DP), Allocatable :: Work(:)
  Integer :: NWork
  Character (len=100) :: foo

  Write(*,*)'Enter matrix size: '
  Read(*,*)N
  M = N

  ALLOCATE(A(N,N), SA(N,N), U(N,N), VT(N,N), Sv(N))
  CALL Random_Number(A)
  A(:,N) = A(:,1)
  Write(*,*)'A = ['
  Do K1 = 1, N
     Write(foo,*)'(1X, ', N,'ES13.5,1A1)'
     Write(*,foo)(A(K1,K2), K2=1, N), ';'
  End Do
  Write(*,*)']'

  CALL SVD(A, Sv, U, VT)
  Write(*,*)Sv

  Write(*,*)'U = ['
  Do K1 = 1, N
     Write(foo,*)'(1X, ', N,'ES13.5,1A1)'
     Write(*,foo)(U(K1,K2), K2=1, N), ';'
  End Do
  Write(*,*)']'

  Write(*,*)'VT = ['
  Do K1 = 1, N
     Write(foo,*)'(1X, ', N,'ES13.5,1A1)'
     Write(*,foo)(VT(K1,K2), K2=1, N), ';'
  End Do
  Write(*,*)']'

  Stop
End Program SVDtest
\end{lstlisting}

\section{Function \texttt{PseudoInverse(A [, Ikeep])}}
\index{PseudoInverse@Function \texttt{PseudoInverse(A [, Ikeep])}}

\subsection{Description}

This function returns the pseudo-inverse of the matrix $A$.

\subsection{Arguments}

\begin{description}
\item[\texttt{A}:]  Real double precision two dimensional array. The
  matrix whose pseudo-inverse we want to know.

\item[\texttt{Ikeep: }] Integer. Optional. Number of singular values
  that we want to keep to perform the pseudo inversion. If not present
  the inverse matrix is computed. 

\end{description}

\subsection{Output}

A two dimensional array. Contains the pseudo-inverse of $A(:,:)$.

\subsection{Example}

\begin{lstlisting}[emph=PseudoInverse,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Using lapack library to compute the inverse.,
                   label=pseudoinverse]
Program SVDtest

  USE NumTypes
  USE LapackAPI

  Real (kind=DP), Allocatable :: A(:,:), B(:,:)
  Character (len=100) :: foo
  

  Write(*,*)'Enter matrix size: '
  Read(*,*)N

  ALLOCATE(A(N,N), B(N,N))
  CALL Random_Number(A)
  Write(*,*)'A = ['
  Do K1 = 1, N
     Write(foo,*)'(1X, ', N,'ES13.5,1A1)'
     Write(*,foo)(A(K1,K2), K2=1, N), ';'
  End Do
  Write(*,*)']'

  B = PseudoInverse(A)

  Write(*,*)'B = ['
  Do K1 = 1, N
     Write(foo,*)'(1X, ', N,'ES13.5,1A1)'
     Write(*,foo)(B(K1,K2), K2=1, N), ';'
  End Do
  Write(*,*)']'

  A = MatMul(A,B)
  Write(*,*)'A = ['
  Do K1 = 1, N
     Write(foo,*)'(1X, ', N,'ES13.5,1A1)'
     Write(*,foo)(A(K1,K2), K2=1, N), ';'
  End Do
  Write(*,*)']'

  Stop
End Program SVDtest
\end{lstlisting}

\section{Function \texttt{expM(M)}}
\index{expM@Function \texttt{expM(M)}}

\subsection{Description}

This function returns the exponential of the hermitian matrix $M$.

\subsection{Arguments}

\begin{description}
\item[\texttt{M}:]  Complex double precision two dimensional array. The
  hermitian matrix to exponentiate.
\end{description}

\subsection{Output}

A two dimensional array. Contains $exp(M)$.

\subsection{Example}

\begin{lstlisting}[emph=expM,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Using lapack library to exponentiate matrices.,
                   label=expm]
Program TestExpM

  USE NumTypes
  USE LapackAPI

  Complex (kind=DPC), Allocatable :: M(:,:), Res(:,:), P(:,:), R2(:,:)
  Real (kind=DP), Allocatable :: V(:)
  Real (kind=DP) :: fac = 1.0_DP
  Integer :: N = 8
  
  Allocate(M(N,N), V(N), P(N,N), Res(N,N), R2(N,N))

  Do I = 1, N
     CALL Random_Number(V)
     M(I,:) = Cmplx(V(:))/10.0_DP

     CALL Random_Number(V)
     M(I,:) = M(I,:) + Cmplx(0.0_DP,V(:))/10.0_DP
  End Do
!  M = Cmplx(0.0_DP)
!  Forall (I=1:N) M(I,I) = Cmplx(V(I))/10.0_DP



  M = M + Transpose(Conjg(M))
  CALL Show(M)
  R2 = expM(M)

  Res = Cmplx(0.0_DP)
  P = Cmplx(0.0_DP)
  Forall (I=1:N) P(I,I) = Cmplx(1.0_DP)
  Do I = 0, 1000
     Res = Res + P/fac
     P = MatMul(P,M)
     fac = fac*Real(I+1,kind=DP)
  End Do
  
  Res = Abs(R2-Res)
  Write(*,*)Sum(Res)
  CALL Show(Res)


  Stop
CONTAINS
  
  Subroutine Show(AA)

    Complex (kind=DPC), Intent (in) :: AA(:,:)

    Write(*,*)'[ '
    Do I = 1, Size(AA,1)
       Write(*,'(1000ES13.5)')(AA(I,J), J=1,Size(AA,2))
    End Do

    Return
  End Subroutine Show


End Program TestExpM
\end{lstlisting}





% Local Variables: 
% mode: latex
% TeX-master: "lib"
% End: 

